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ABSTRACT 


An interim report is given of a calculation aimed at re- 
constructing the three-dimensional distribution of optically 
thin EUV emission sources associated with solar active regions, 
from two-dimensional observations (projections) recorded by 
the Goddard spectroheliograph on the OSO-7 satellite. Its 
relation to other image reconstruction methods in the literature 
is briefly discussed, together with special requirements im- 
posed in the solar case such as a knowledge of the true solar 
rotation function. A useful correlation criterion for es- 
tablishing the physical validity of solutions is given. 
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INTRODUCTION; 


The Goddard X-ray and extreme ultraviolet spectrograph on the seventh 
Orbiting Solar Observatory scanned the sun and Its corona during the 
period October 2, 1971 to July 1974 » producing spectrohellograms In 
selected wavelengths corresponding to known, optically thin, emission 
lines emitted by coronal Ions [1]. In addition to "large raster" spec- 
trohellograms covering the whole sun, "small raster" spectrohellograms 
about five arc-minutes on a side were obtained which provided intensity 
distributions for particular active regions. A typical small raster con- 
tour plot Is shown in Figure 1. In the case of many active regions, a 
large number of such rasters were obtained at a given wavelength over a 
period of days as the region rotated across the solar disk. The avail- 
ability of these successive images or projections prompted the present 
study to determine whether the three-dimensional emission distribution 
within an active region can be reconstructed from a pair of OSO-7 rasters. 
This article Is an Interim report on the calculation, which during the 
past year has been carried to the point where It may be possible to de- 
termine the feasibility, or otherwise, of such a reconstruction from the 
available data. The calculation was developed Independently of similar 
methods in other fields (e.g. the medical field, in which striking ad- 
vances have recently been reported (reference [2], and other references 
in the same issue of IEEE Trans, on Nuclear Science)), but relationships 
will be discussed below. Two difficulties are encountered here beyond 
those dealt with by Image reconstruction techniques in the laboratory. 
First, the source distribution is not necessarily constant with time, 
as well as possessing no distinct boundaries either within or without. 
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Figure 1, The small square projected here on the solar disc represents a typical 
0SO7 small raster spectroheliogram observed in the Fe XV emission 
line wavelength 284. lA. This particular raster, obtained on January 
22, 1972, shows intensity contours of the active region McMath No. 

116< 


Second) the solar rotation rate is not well known at coronal levels. 

Thus the general difficulty of the problem has been greater than anti- 
cipated. In the course of the calculation) however, a three-dimensional 
lattice of points is constructed which can be superlo 4 >osed on an active 
region at different times and has proved to be useful for the Investi- 
gation of its rotation. 

Much of the calculation consists of basic geometrical and trans- 
formational steps which will be included for completeness. 

GIVEN DATA; 

The original OSO-7 measurements consist of a 23 x 16 Irregular 
array of Intensity elements over a total field of view of approximately 
5x5 arc minutes for each small raster, which requires about one minute 
of time to record. For non-flaring active regions, up to twenty indivi- 
dual rasters are melded together to improve statistical accuracy. Then 
the elements are Interpolated to produce a uniform array of 2500 Intensity 
values 1(1, j), 1 and j each taking values from 1 to 50, with an equal 
angular spacing of 6 arc seconds between points. Two such melded rasters, 
recorded at successive times such that a suitable solar angle separates 
them, are used In the following analysis and will be referred to as 
rasters A and B. Associated pertinent Information is included with each 
of these melded rasters on a computer tape, such as: heliocentric co- 

cordlnates of the central point (r^, 0^) and four corner points (rj^, 9^^; 
^2* ^2* ^3* ^4> ^ 4 ) * times of beginning and end of the period covered 

by the Individual rasters Included In the meld; solar pole position 
angle P and sub-earth solar longitude Bq. 
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PROCEDURE: 


A basic coordinate system, used throughout the calculation, was 
constructed as shown In Figure 2 , with the sun's center as origin, x- 
axls passing through the earth, and z-axls directed toward an arbitrary 
direction here taken to be the projection of celestial north. (The 
computations are somewhat simplified If the projection of the solar 
north direction is taken as the z-axis, since P is then zero throughout, 
but the calculation is presented here as it was carried out.) The 
rasters A and B are considered to be two-dimensional images lying in the 
(y,z) plane. First raster A will be dealt with. 

1 . Use of raster A; 

The coordinates of the center point and corner points of raster 
A are given by: 

(Xc* Zc^ * sin( 0 j. + P), rj. cos( 0 ^ + P)) 

(yi» zi) “ (-^1 sin( 0 i + P), cosCGj + P)) 

(X 2 * Z 2 > “ <-*^2 sin(e^ + P) , V 2 cos (92 + P)) ( 1 ) 

(X3, Z3) - (-r3 sin(03 + P) , r3 cos (©3 + P)) 

(X4> 24) » (-r4 sin(04 + P) , cos(04 + P)) 

Using these coordinates a two-dimensional grid of x N2 points (yp Zp q) 
lying within the raster is constructed in the following way. The lengths 
of the sides of the raster are given by - [iy2 ' 3^3) (*2 ” * 3^^1 » 

^2 “ t^X4 " X3)^ (24 - 23)^1^ while the slopes of the sides are 

raj^ “ (Z3 - Z2)/(y3 - y2> and m2 “ (23 - - y^) . Choosing the 

distances between the grid points as A-f.3 = 1), AI2 “ ^2^2^(^2 ^ 

with K], and K2 as desired fractions of the raster dimensions to be covered 
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Figure 2. Coordinate system used in the text, with origin at the center of the 
Sun, z axis in the direction of celestial north, earth observer at the 
point (d, O, O) on the x axis. D represents the solar rotation axis, 
of which the projection on the (y, z) plane lies at the angle P to the 
north direction; ^ is tilted out of the (y, z) plane by the angle Bq, A 
detailed discussion of the solar parameters P and Bq may be found 
on pages 307-9 of reference (8). R represents an observed point in 
an EUV raster, lying in the (y, z) plane and defined by the solar disc 
coordinates (p,^). 
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by the grid, a £irst grid point is constructed as: 
yi 1 - sign(y3 - y4>[^(N2 - 1) ^£ 2/(1 + 

+ sign(y 3 - y 2 )(^(Ni - 1 ) ( 2 a) 

^ 1,1 " (similarly) 

Points (yp 3 ;zp^l) are then constructed as 

Yn 1 “ sign(y 2 - y 3 )[(p - 1 ) ^£1/(1 + m^"^)^] + 1 

P,i (2b) 

Zp 3 “ (similarly) 

and general points (yp^q'.Zp^q) are further constructed according to 
yp,q sign(y 4 - y 3 > t(q - 1 ) + m 2 "^)^] + yp^l 

( 2 c) 

Zp^q “ (similarly) 

This produces the desired two-dimensional grid, centered within raster A. 
Corresponding polar coordinates of these points in the (y,z) plane are 


Pp,q ■ 'yp.q^ 

0p_q - tan-1 [-yp^q/zp^ql 


( 3 ) 


These grid polar coordinates are next used to produce a three-dimensional 
lattice of x N 2 x N 3 points (Xpqj-, ypq^.> Zpqj.) extending throughout 
the volume of solar atmosphere viewed by the raster, consisting of 
successive two-dimensional grids separated by a chosen distance li 


dp2p_q + d[c^ „ (d2 + p2„ „)((! + rZ)2 - „)) 


p,q 


p.q 


p»q 


p.q,r 


d2 + o2 


p,q 


S-p.q.r - -(tp.qW - Xp.q.r^ *1" »p,q'''‘ 
'p.q.r * 'dp.q - Xp.q.P I'p.ql'"* 


( 4 ) 


d being the earth-sun distance in solar radii. 
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In this lattice, successive points (p.*i,r), (p,q> r + 1) , (Ptq, r 2), ... 
lie along a line of sight joining the earth (observer) to the point 
(yp^q; lying in the (y.z) plane (see Figure 3); thus the lattice, 

conposed of x N2 such lines converging to the earth, is not exactly 
a cubic lattice but departs from parallelism by the small angle of less 
than 1 arc*minute or 0.0003 radian 'existing between the lines. For the 
present purpose this small departure can be neglected. 

The solar latitude Bpq^ and the hei^t Hpq^ above the surface of each 
of the lattice points is given by 

®P,q,r ■ (1 + ®i‘'"M’«p,q,r ®o - yp,q,r ^ ®o 

( 5 ) 

+ Zp^q^j. cos P cos Bq] 

Hp^q^j. • (6.9598 x 10^)rjf (in km) 

The volume of solar atmosphere enconq>assed by this lattice may be 
considered to be composed of discrete rectangular source blocks, each 
centered on a lattice point (Xpq^, Fpqr» ‘pqr^* that along a line of 
si^t (p,q) the total intensity Ip,q emitted toward the observer is the 
sum of separate contributions S(p,q,r), S(p,q, r 1), S(p,q, r + 2) ... 
emitted by the source blocks centered on the corresponding lattice points , 
thus Ip^q emerges from the end of a source column of rectangular cross- 
section whose axis is the line of sight (p,q). To obtain Ip^q from the 
actual raster intensities I(i,j), the latter values are summed over values 
of i,j which lie within the cross-section of the source coluom. In this 
step the following equivalence between raster indices (i,j) and rectangular 
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Figure 3. This figure illustrates, in a schematic way, the relation of the lattice points given by eouations 
to the original grid points (2). 



coordinates (y,z) is needed: 


y<i.j)=-^ty3+^ C2J-l)(y4-y3)l+a--^)[y2+^ (Zj-lXyi-yj)] 

z(i,j) = (similarly) 


( 6 ) 


Now, by transforming from the triplet of variables p,q,r :o two new 
variables t,i by the transformation 


t = N 2 <p - 1) + q 
i = N 2 N^(p - 1) + N 3 (q- 1 ) + r 


(7) 


a set of linear equations is constructed which expresses the fact that 

each intensity value L - = lx corresponds to a source column t composed 

P ? 4 ^ 

of contributions S^^i, 


1(1) = S(l,l) + S(l,2) + ... + S(1,N3) 

1(2) = S(2,N3 + 1) + S(2,N3+2)... + S(2,2N3) 


( 8 ) 


I(N^ N2) = S(N 3 ^ N2, \ ^2 % ■ ^3 + 1 )+. ..+S(N;L ^2 N 3 ) 

2. Use of raster B; 

Since the active region registered on raster A at time t^ has 
rotated through an angle 6 ^ 360 ((tg - t^)/24) about the solar axis by the 
time it is registered on raster B at time tg, one must rotate the lattice 
blocks centered on (Xpqj., ypqj-> ^pqj.) by the same amount and around the 
same axis . The rotation of a vector r through an angle 6 about an arbitrary 
axis D, in a clockwise direction when viewed along the positive direction 
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of D, is given by the transforoiation [3]: 

r' - r + sin 6 (D x r) + (1 - cos 5)(D x (D x t)) (9) 

Application of this rotation to the original lattice points 

Ypqj., Zpqr) results in the new lattice (Xp^^., Ypqj.. ^pqj.) superimposed 

on the active region in its later position, as registered in raster B: 

^,q,r *pqr ^ft)‘Sqr ‘ VD*ypqr) ^pqr 

+ (1-cos 6pqj.) [^(op.ypqr ■ ^*Xpqr) ~ YD^Yo’^pqr " °^D*^pqr^^ 

^P>q»r “ ypqr ^YD'^pqr ’ “'D’^pqr) ^pqr 

»• (1 - cos 6pqr^ tY|)(ft)*^pqr ~ YD*Ypqr^ ” ®D^®D*ypqr ” A)*^pqr^^ 

^P>q»r * ^pqr (®D*ypqr “ ft)*^pqr^ ^pqr 

+ (1-cos 6pqr)[ofD(YD.Xpqi- - OD-Zpqi.) - A)^^‘^pqr ' M}*ypqr^^ 

( 10 ) 

where = sin Bq, ft) = - sin P cos B^, = cos P cos and the rotation 

angle 6pqj. - a function of the latitude Bpqj. of each point - has been 
taken here initially to be a function of the form a+ bsin Bpq^, based 
on the latitude dependence of solar rotation which is generally accepted 
in the literature. Remarks on the importance, for this calculation, of 
finding a correct function to represent 6,jqj. will be made in concluding 
sections . 

The corners of each source block centered on a given lattice point 
(Xpq^, ypqr» ^pqr^ similarly rotated into their new positions. With 


10 


the aid of the new corner coordinates, the equations of the six bounding 
planes of a given source block may be obtained; for convenience we 
denote them here as : 

aix + bj^y + c^z + d]^ =0; a4x + b^y + c^z + d^, = 0 

32X + h2y + C2Z + d2 = 0; ajx + b^y + c^z + d^ = 0 (11) 

33X + b3y + C3Z + d3 = 0; a^x + b^y + c^z + dg = 0 

Now, unlike the simpler situation existing in the case of raster A, 

a line of sight to the point (Xpqj., Ypqj., Zpq^ does not intersect any 
other point (Xp'q'j.'* ^p'q'r^’ ^P^q^r^^ general, so that a corresponding 
observing cylinder of rectangular cross-section will intercept the source 
blocks along its path in various fractional volumes. Thus such an ob- 
serving cylinder (p,q,r), corresponding to an observed intensity Ip q j., 

will intersect a certain volume V_ _ _ of the source block centered on 

P ^ 

(Xpqr» ^pqr» ^pqr^ » intersect other source blocks (Xyvw>'^uvw»2uvw) 

in the vicinity of the line of sight in further particular volutins 
If the full volume of a source block is V^, then the contribution of 
source block (Xpqj., Ypqj-, Zpqj.) to the observed intensity Ip^q^r 
Vpqr/Vjj, and similarly the contributions of neighboring blocks (u,v,w) 
are V^^/Vq. These contributions will then all be less than unity. 

The four planes bounding a typical observing cylinder may be written 
for conciseness as; 

ayx + b^y + Cyz + dy “0; aqx + bgy + cqz + dg = 0 

( 12 ) 

agx + hgy + cgz + dg = 0; a^QX + bj^Qy + c^qz dj^Q = 0 
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To obtain the actual value of the volume of Intersection of observing 
cylinder and source block - specified by the intersection of the set of 
planes (11) with the set of planes (12) - a Monte Carlo procedure is used 
in which random points are thrown into the source block, and the numbers 
of these points which also fall within the observing cylinder are counted. 
(This turns out to be a difficult problem to solve analytically.) 

The contributions Vpqj./V^ and are just the coefficients 

^t i 3 linear equations which state that the observed 

intensities (obtained by appropriately summing intensities I(i,j) 

of raster B over the cross-section of observing column (p,q,r), as was 
done for raster A) are sums of source contributions along the line of 
sight to the point (p,q,r), and which will be added to the first set con- 
structed for raster A. To formulate this additional set we transform from 
the variables (p,q,r) to the variable t by the relation: 

t = Nj N2 + N2 N3 (p- 1) + N3(q-1) + r (13) 

with the relation between i and (p,q,r) remaining as in (7). This results 
in the set : 

I(N]^ N2 + 1) = ^ a(N^ N2 + l,i) S(N^ N2 + l,i) 

I(N^ N2 + 2) = Z a(N^ N2 + 2,i) S(Ni N2 + 2,i) (14) 

I(N^ N2 + N2 N3) = Z a (N^ N2 + N2 N3, i) S(Ni N2 + % N2 N3,i) 

with i taking values from 1 to N2 N2 as before. 
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These equations, together with the first set ( 8 ), can be represented 
by the following full matrix equation for the unknown source vector S 
whose elements are the emission strengths of the individual source blocks; 

i - |Als , 05) 

in which the coefficient matrix |a|, a sparse matrix, is of dimensions 
(Nf N 2 + Nj^ N 2 N 3 ) X (Nj^ N 2 N^). The general form of |a| can be repre- 
sented schematically as in Figure 3, 

DISCUSSIC»J AND PRELIMINARY RESULTS: 

1. Relationship to other three-dimensional reconstruction methods: 

Gordon and Herman [4] classify existing three-dimensional reconstruction 
methods into four main categories; sunmiation, use of Fourier transforms, 
analytic solution of integral equations, and series expansion approaches. 
Budinger and Gullberg [2] expand this classification further into thirteen 
detailed procedures, one of which consists of direct matrix and linear 
equation methods (included within the series expansion approaches described 
by Gordon and Herman). The present calculation falls into this category 
which is discussed in a physically meaningful way by Budinger and Gullberg. 
Further, it deals with a "situation where we are limited in the number 
of views" in which Budinger and Gullberg note that the resulting matrix 
may be close to singular, and may therefore require for its solution the 
use of iterative methods. 

Direct inversions of the coefficient matrix |a| have been attempted, 
as a first trial, for pairs of rasters separated by solar rotation angles 
from 15® to 47®. The three-dimensional lattice used was 5 x 5 x 5, giving 
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General appearance of the final coefficient matrix A; the non-zero 
elements denoted by x have values between zero and unity. All 
elements other than unity or x are zero. 
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dimensions of 125 x 150 for |a| . Unique solution vectors S, in the 
mathematical sense of satisfying the original equation (15), are obtained 
which however possess some negative - nonphysical - elements. These 
solutions may nevertheless contain some physical reality. (In this 
connection one may note that a standard procedure in some of the ART 
(algebraic reconstruction techniques) algorithms [2], [4] is to set 
negative elements equal to zero.) We are now at a stage in the calculation 
in which it may be possible to assess the reality, or non-reality, of 
these solutions. A correlation method which we expect to use for this 
purpose, and whose general principle does not appear to have been 
suggested previously, is described in the following section. 


2. A Correlation Criterion for Validity of the Solution: 

When a solution S(xpq^, ypqj- » ^pqr^ matrix equation (15) is 

obtained, one would like to establish whether it is a physically valid 
solution or not. This is not a trival problem in the case at hand, since 
active region distributions are expected to be diffuse rather than compact, 
as mentioned in the Introduction where it was noted that in this respect 
the present problem is more difficult than the medical three-dimensional 
reconstruction problem. A correlation method is described here which 
may be of general interest in such cases. The method consists of running 


the calculation, on a given pair of rasters, with successively decreasing 
grid sizes which produce corresponding solution sets S^(x^pq^, y^pq^.. 


zW>> ^‘’pqr^ • ^ Slven 


pair of solution sets, say and , one first interpolates the function 


S®(x®__^, y®___, z^nnT-) find its values S^(x^__„, y^ 
pqr’ pqr’ pqr' ' pqr* 


pqr' pqr^ 
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the lattice points of the second solut:.on For a given r - i.e. for 

a given height in the solar atmosphere - the following normalized 
correlation coefficient is then constructed: 


r^(a,b) 


Ptq 




S^Cx* 


pqr» y pqr’ 


z pqr^ 


(16) 


2^q[S><x^p^j., y^pq^, *^pqr^^ 

where represents the larger of or S^. This coefficient is seen to 

approach unity the more closely the solution field resembles the 

solution field The final overall parameter describing the correlation 

of solution with is then r(a,b) = T. r^(a,b). The hypothesis under- 

r 

lying the use of this parameter is that a physically valid solution field 
will be Independent of the size of the lattice used. It is expected to 
be of value for example in improving the trial rotation period used in 
the calculation, since one can vary the assumed rotation period so as to 
maximize r(a,b) or more generally the quantity - S r(i,j). 


3. Requirements imposed by the data: 

The present procedure could be extended to include more than two 
projections (rasters) without undue difficulty, but two projections only 
were used in this first attempt because of limitations imposed by the nature 
of the solar sources dealt with. An ideal source region, for the purpose 
of this inversion calculation, would possess the following properties: 

(a) the source distribution should remain reasonably constant with time; 

(b) its emission should be much greater than that of any background 
which is present; (c) its size should be appreciably less than the 
raster size, i.e. than the cross-section of the observing "beam"; and 
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(d) its rotation rate should be accurately known. Actual solar active 
regions fall short of the ideal in all of these respects. The solar 
rotation rate is known only to a first approximation, being a function 
of several variables such as solar latitude, altitude above the solar 
surface, type of feature considered (sunspots, coronal features, magnetic 
features, etc.). Recent papers which give further references on this 
topic are those of Henze and Dupree [5] and Simon and Noyes [6]. The 
present uncertainty which exists in the field may be illustrated by the 
fact that while Simon and Noyes find "constancy of rotational rate be- 
tween chromosphere and corona" and Henze and Dupree similarly "cannot 
discern a variation in rotation rate between chromosphere and corona", 
El-Raey and Scherrer [7] find that "the upper chromosphere and lower 
corona rotate from 5 to 8Z faster than either the photosphere or corona". 
On the other hand, since the present calculation requires an accurate 
solar rotation function as input for a successful solution, it also has 
the potential of establishing such a function unambiguously. 

The further limitations imposed by conditions (a), (b) and (c) above 
are presently being explored in the OSO-7 data. 
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